Eddington-inspired Born-Infeld gravity. 
Phenomenology of non-linear gravity-matter coupling 
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Viable corrections to the matter sector of Poisson's equation may result in qualitatively different 
astrophysical phenomenology, for example the gravitational collapse and the properties of compact 
objects can change drastically. We discuss a class of modified non-relativistic theories and focus on 
a relativistic completion, Eddington-inspired Born-Infeld gravity. This recently proposed theory is 
equivalent to General Relativity in vacuum, but its non-trivial coupling to matter prevents singular- 
ities in early cosmology and in the non-relativistic collapse of non-interacting particles. We extend 
our previous analysis, discussing further developments. We present a full numerical study of spher- 
ically symmetric non-relativistic gravitational collapse of dust. For any positive coupling, the final 
state of the collapse is a regular pressureless star rather than a singularity. We also argue that there 
is no Chandrasekhar limit for the mass of non-relativistic white dwarf in this theory. Finally, we 
extend our previous results in the fully relativistic theory by constructing static and slowly rotating 
compact stars governed by nuclear-physics inspired equations of state. In the relativistic theory, 
there exists an upper bound on the mass of compact objects, suggesting that black holes can still 
be formed in the relativistic collapse. 

PACS numbers: 04.50.-h, 98.80.-k 



I. INTRODUCTION 

The beauty and the beast of Einstein's General Rela- 
tivity (GR) are encoded in the non-linearity of its field 
equations. Already in vacuum, GR describes the dynam- 
ics of non-linear objects, like black holes [l|. In order to 
study the formation of black holes in dynamical situa- 
tions, e.g. during a stellar collapse, one needs to couple 
the vacuum theory to matter. How to include this cou- 
pling in a proper way was one of Einstein's main con- 
cerns. The requirement of stress-energy tensor conserva- 
tion, V p T' il/ = 0, which in turn implies geodesies motion, 
together with Bianchi's identities, V = 0, naturally 
suggests a linear coupling between the Einstein tensor 
and the stress-energy tensor, G^ v oc T^ v . Indeed, under 
quite generic assumptions (see e.g. 0) Einstein's equa- 
tions are the most general field equations which involve 
the stress-energy tensor linearly. Nonetheless, it comes 
as a surprise that a highly non-linear theory as GR is 
just linearly coupled to matter. In this paper we inves- 
tigate modified theories of gravity which account for a 
non-linear coupling to matter, while retaining many ap- 
pealing features of Einstein's theory. 

In the weak-field regime, GR has brilliantly passed all 
experimental and observational scrutinies so far [H, Q], 
but new effects are still viable in the strong-field, non- 
linear regime. Furthermore, while current experiments 
have confirmed the weak-field behavior of GR in vacuum 
or in orbital motion, performing null tests of Einstein's 
theory inside matter may be extremely challenging be- 
cause, due to the equivalence principle, purely gravita- 
tional effects are hard to disentangle from those due to 



non-standard matter coupling. In fact, one usually as- 
sumes a minimal coupling, then solves the full Einstein 
equations in some relevant situation - e.g. in cosmolog- 
ical settings, in the interior of Sun-like stars or inside 
more compact objects like neutron stars (NSs) - and fi- 
nally compares theoretical models with observations. 

On the other hand, GR suffers from severe theoretical 
problems and long-standing observational puzzles, which 
may be precisely related to our poor understanding of the 
gravity-matter coupling. For example, the dark matter 
problem (see e.g. Ref. |5| for a review) may be explained 
by invoking new fundamental interactions [f| 0] , rather 
than assuming exotic particles. Similarly, the cosmologi- 
cal acceleration of the universe may be explained in terms 
of more complicated interactions, rather than postulat- 
ing the existence of a mysterious form of dark energy (see 
e.g. Ref. for a recent review on cosmology in modi- 
fied gravities). These postulates are somehow the modern 
versions of the aether, suggesting that perhaps something 
is missing in our understanding of gravitational interac- 
tions inside matter. 

Likewise, it is well known that the dynamical evolu- 
tion of matter fields in GR is generically plagued by the 
formation of singularities (e.g. the Big Bang or those 
forming in the gravitational collapse of matter fields), 
which signal a break-down of the theory. According to 
Penrose's cosmic censorship Q, these singularities must 
be covered by an event horizon, i.e. a black hole must 
form as an outcome of the gravitational collapse. How- 
ever, the cosmic censorship remains a conjecture and its 
validity (or even its precise formulation) is an open issue. 

In a spherically symmetric stellar collapse, locally 
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naked singularities can be generically formed [1CH12I ] . Ini- 
tial density gradients can produce shear in the collapsing 
fluid which, in turn, delays the formation of an appar- 
ent horizon, leaving the singularity locally naked [ 1 3| . In 
some cases, the singularity can even be globally naked, 
i.e. connected to distant observers by timelike or null- 
like geodesies, thus s ugg esting that the censorship can 
be evaded (see Refs. [14l - [l6| for detailed discussions on 
this topic). However, these "shell- focusing" singularities 
may be limited to the spherically symmetric case (l7j . 
Relaxing the assumption of spherical symmetry, the sit- 
uation is even more controversial. While Shapiro and 
Teukolsky reported numerical evidence for formation of 
naked sin gula rities in the collapse of highly prolated gas 
spheroids [18| , Wald and Iyer [l9[ pointed out that similar 
properties - namely the absence of trapped surfaces ly- 
ing on their maximal slices in a portion of singular space- 
time - are also present in a Schwarzschild spacetime with 
a particular choice of the time slice. These studies show 
that, although there appears to be growing evidence in 
support of the cosmic censorship, its validity in GR re- 
mains an open issue which is far to be solved [13, H3| ■ 

However, we stress here that the issue of singularity 
formation and the related cosmic censorship strongly de- 
pend on how gravity interacts with matter. Thus, the 
outcome of a gravitational collapse may be different if 
the gravity-matter coupling is modified. 

In this paper, following previous works [2l| - [26| . we take 
a different perspective and investigate phenomenologi- 
cally viable modifications to GR that account for non- 
trivial coupling to matter while being equivalent to Ein- 
stein's theory in vacuum. As a consequence, this class 
of corrections results in a qualitatively different phe- 
nomenology once matter is included (see also Ref. |27|). 
For concrcteness, we will focus on a recent proposal by 
Banados and Ferrei ra [25l| (see also some previous at- 
tempt by Vollick [2l|-[23|Twhich we shall call Eddington- 
inspired Born-Infcld (EiBI) gravity, not to be confused 
with Eddington-Born-Infeld [28| or Born-Infcld-Einstein 
theory [29| . While being conceptually different from each 
other, these theories share the same determinantal form 
of Born-Infeld non-linear electrodynamics [3(3, l3l| - In 
this sense, similarly to non-linear electrodynamics, EiBI 
gravity can be thought as an effective prototypical theory 
in which resummation of higher order curvature terms 
- which are qualitatively similar to those expected by 
a quantum gravity completion - is advocated in order 
to resolve curvature singularities [32| . It turns out that 
these corrections effectively account for some non-linear 
matter coupling to usual Einstein's gravity. 

Interestingly enough, such corrections can even affect 
the non-relativistic gravitational regime. We shall discuss 
the non-relativistic phenomenology in detail and most of 
our results apply to any theory whose non-relativistic 
limit is described by the following modified Poisson equa- 
tion m 



V 2 $ = 4irGp + -V 2 p, 



(1) 



where $ is the usual gravitational potential, G is the 
gravitational constant and p is the matter density. 
Within the parametrized post-Poissonian approach pro- 
posed in Ref. [33j , the equation above is the most general 
spatially covariant Poisson equation, which is first order 
in $ and p and reduces identically to standard Laplace's 
equation, V 2 $ = 0, in vacuum. While possible second or- 
der terms are strongly constrained by tests of the equiva- 
lence between gravitational and inertial mass to one part 
in I0 12 (see e.g. Ref. [HI and references therein), the term 
proportional to n in Eq. (TTJ) is presently mild constrained. 

To our knowledge, null tests of Poisson equation inside 
matter, in order to constrain the extra term in Eq. ([1}, 
have not been conceived yet. Mild observational con- 
straints on k come from solar physics observations [33j . 
\k\ < 3 x 10 5 m 5 s 2 /kg. This bound has been derived 
by comparing the observed solar neutrino fluxes and he- 
lioseismology observables to the predictions of modified 
solar models. Such constraints are quite mild, due to 
some inherent uncertainties on the Sun interior. Very 
recently, considerabl y st ronger observational constraints 
were derived in Ref. [35[ , provided direct measures of the 
central density or of the radius of a NS are available (see 
also Ref. [26|). However, some uncertainty exists on the 
interior of compact stars and usually GR is assumed to 
infer the magnitude of the central density, given some 
"observable" quantity such as the mass and the radius. 
Thus, table experiments in some controlled setting would 
be highly desirable. 

The scope of this paper is twofold. First, we wish to 
discuss the phenomenology of Eq. (JTJ) in relation with 
stellar collapse and static stellar configurations in the 
non-relativistic limit. Our purpose is to show that, while 
being compatible with current observations, a simple 
modified Poisson model leads to a new interesting phe- 
nomenology, even for arbitrarily small values of k. Sec- 
ondly, we shall discuss a well-motivated relativistic com- 
pletion of Eq. (P), EiBI gravity [n|. In Ref. [HI, we 
reported some results on the gravitational collapse and 
compact star solutions in this theory. Here, we give fur- 
ther details and discuss new developments. 

The plan of this work is the following. In Sec. |TT] we 
briefly review some known and new features of EiBI grav- 
ity [231 . Sections |HT] and [IV] are devoted to the phe- 
nomenology of the modified Poisson equation ([1} , which 
includes the non-relativistic limit of EiBI gravity. Sec- 
tion |HT] is devoted to solve the hydrodynamics equations 
describing the non-relativistic collapse of non-interacting 
particles. In pass, we shall correct a wrong result in 
Ref. [1^ which, however, docs not affect the final pic- 
ture. In Sec. II VI we show that the end-point of the grav- 
itational dust collapse is a regular, pressureless star and 
we discuss other non-relativistic stellar equilibrium con- 
figurations. Finally, in Sec. |V]) we discuss relativistic 
stellar models using nuclear-physics based equations of 
state. In Sec. IVII we draw our conclusion and discuss 
open issues and future developments. 
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II. EDDINGTON-INSPIRED BORN-INFELD 
GRAVITY 

EiBI gravity is described by the following action [25| 



K 



d 4 x 



g ab + nR ab 



Sm [g, # 



(2) 

where Sm [g, ^m] is the matter action, ^m generically 
denotes any matter field, R ab is the symmetric part of the 
Ricci tensor built from the connection T^ b and A is related 
to the cosmological constant, A = (A — 1)/k, so that 
asymptotically flat solutions are obtained when A = 1. 

EiBI gravity is reminiscent of a Born-Infeld action for 
non-linear electrodynamics (30L l3lj , where here the Ricci 
tensor plays the role of the field strength F ab . Moreover, 
the metric g and the connection V are considered as in- 
dependent fields [25|, as in Palatini's approach to GR. 
Similarly to f(R) gravities (see e.g. [36|) the Palatini for- 
mulation is not equivalent to the metric one. In the met- 
ric approach, EiBI theory contains ghosts, which must be 
eliminated by adding extra terms to the action [U [2!| . 
Furthermore, it is assumed that matter minimally cou- 
ples to the metric tensor only, i.e. the matter action Sm 
in Eq. ([2]) only depends on the metric g and on the matter 
fields, but not on the independent connection T. 

The field equations are conveniently written as [25| 



I 



r ab(<?) = ^<f d {d a qbd + d b q ad - d d q ab ) , 




KR ab (r) , 
b _ KT ab\ 



(3) 

(4) 
(5) 



where q ab is an auxiliary metric, T ab is the standard 
stress-energy tensor and i? a f,(T) is the symmetric part 
of the Ricci curvature tensor. Crucially, since matter is 
minimally coupled to the metric g ab , the field equations 
above imply the usual conservation of the stress-energy 
tensor, \7 a T ab = 0, where V is the covariant derivative 
written in terms of the metric tensor g. Note that q ab is 
the matrix inverse of q ab and differs inside matter from 

q cd g ac gbd- 

When T ab = 0, the auxiliary metric q coincides with 
the physical metric g. Thus, in vacuum T^ c is the metric 
connection. Indeed, in absence of matter fields EiBI the- 
ory is equivalent to GR [25j . In presence of matter, the 
small n limit of the field equations reads 



R a b(T) — T ab — -Tg ab + K 



S a b ~^Sg ab 



+ 0(k 2 ). (6) 



where S ab = T c a T cb ~^n ab . 
ory is recovered as k — > 0. In the field equation above, 
two qualitatively different corrections to GR are manifest 
at O(k): those depending on S ab , which are quadratic in 
the matter fields, and those hidden in R ab (T), which im- 
plicitly depend on derivatives of matter fields. The lat- 
ter survive even in the non-relativist limit of the theory, 
which is described by a modified Poisson equation ([1]) . 



\TT nb . Notice that Einstein's thc- 



Since the auxiliary metric q ab is algebraically related to 
the physical metric g ab , EiBI gravity does not introduce 
any extra dynamical degree of freedom with respect to 
GR. In this sense, the theory is similar in spirit to f(R) 
gravities in the Palatini approach [36j . 



A. Two formulations of EiBI gravity 

Let us now investigate the field equations of EiBI the- 
ory in more detail. Written in terms of the new rank- 
two auxiliary tensor q, the theory is reminiscent of a 
particular bi-metric theory (see e.g. Refs. [13, HI] for 
some specific attempt) which, in vacuum, degenerates 
into a single metric theory, GR. On the other hand, in- 
side matter the two metrics are different and the connec- 
tion is metric with respect to q, and not g, cf. Eq. (|3|). 
Nonetheless, the conservation of the stress-energy tensor 
- and geodesic motion therein - is guarantee by a min- 
imal gravity-matter coupling in the action ([2]). This is 
precisely the essence of the theory. 

It should be noted however that there are at least two 
ways of interpreting the theory to which we shall refer to, 
with some abuse, as the "Jordan frame" and the "Ein- 
stein frame" , in accordance to scalar-tensor theories. In 
the Jordan frame matter is minimally coupled to the met- 
ric g and Einstein equations are modified according to 
Eq. ([6]). In the Einstein frame the metric q is treated as 
fundamental and matter is non- minimally coupled to it. 
This appears quite naturally in the small n limit of the 
theory, Eq. ^ , where g^ can be written solely in terms 
of q^y and its derivatives by using Eq. ((4]) . 

Let us however start by pointing a particular symmetry 
of the action. In the Jordan frame, the action reads as 
in Eq. @. In what we define as the Einstein frame, the 
action can be written explicitly by defining a new metric 

lab = g a b + nR ab (r), 



S-y 



-2A 



4 X (^J\lab - nR a b(r)\ - ^n/-^ 
SMhab-KRab(r),VM} • (7) 



Written in this form, the metric 7 is now non-minimally 
coupled to matter. Furthermore, the field equations im- 
pose 7 a b = q ab , but the connection T is not anymore 
the g-mctric connection, since T enters explicitly in the 
matter action through i? a t>(r). 

Neglecting the matter action in S g and 5 7 , the follow- 
ing duality holds 



g a b <->• 7ab . 



A-B-i/A, r^r. 



(8) 



up to an overall factor A. This duality can be seen in 
the vacuum field equations, where the rescaling follows 
from the algebraic relation (|5|). Note that only in the 
case A = 1, i.e. for a vanishing cosmological constant, 
lab = gab- This also shows that the case A = I is spe- 
cial not only because it corresponds to asymptotically 
flat space in both frames, but it is a fixed point of the 
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duality (|5J). If A = 1, negative values of k in the Einstein 
frame correspond to positive values of k in the Jordan 
frame, provided one modifies the matter coupling. 

However, the coupling to matter breaks the symmetry 
above, introducing a non-minimal coupling in the Ein- 
stein frame. Interestingly, the action ([7]) in the Einstein 
frame reduces, in the small k limit, to Einstein's theory 
with non- minimally coupled matter (see also Ref. [27j). 



A / cFxy/^i (R - 2 



.A- 1 



1 



+ / <?Xy/=j 1 - -KR(T) C m {lab - KR ab {Y)) 



(9) 



where we have defined R = j ab R ab - 

Furthermore, if the matter Lagrangian is constructed 
in terms of the metric g, the second term in Eq. ^ in- 
volves other k corrections and some higher order coupling 
to the Ricci tensor and to the independent connection. 
For example, let us consider the small k limit of a mini- 
mally coupled free scalar field in the Jordan frame 



C m {.g a b, 4>) = g a bd a 4>d b (j) . 



(10) 



In the Einstein frame, the matter Lagrangian reads 

£ m (7a6,r,0) = la bd a <j>d b <j>-KR ab {T)d a <pd b <j>, (11) 

leading to the following matter action in the Einstein 
frame 



Sn 



d x \/— 7 



(12) 



-K[R ab + -R lab ) d^d^ + Ci^), 



Note that the action above involves second derivatives of 
the scalar field in the gravity equations hidden inside the 
Ricci tensor evaluated on shell. 

A similar transformation as the one presented above 
persists in the non-relativistic limit. Indeed, if we define 
<f> = $ + jp, Eq. (JTJ) reduces to the standard Poisson 
equation, 



4nGp. 



(13) 



However, in this case the gravitational force would read 
F = V$ = — j V/0, being thus dependent on the 
matter field. 

The modified Poisson equation (fTJ) has an interesting 
interpretation in terms of standard Poisson equation with 
a modified source term. The second term on the right 
hand side of Eq. (fTJ) leads to an effective force: 



F h 
P 



-4 V '' 



(14) 



from which we obtain that P e ff = np 2 /8. Therefore, the 
modified Poisson equation is formally equivalent to the 



standard one supplemented by an effective polytropic 
fluid with equation of state (EOS) P(p) = Kp 1+1 / n , 
where the polytropic index n = 1 and K = k/8. In fact, 
a similar property exists also in the relativistic case. In 
the weak field limit of EiBI theory, one can show that 
an effective pressure and effective internal energy ap- 
pears. This will be explained in detail in a forthcoming 
work SI. 



B. Linear structure of EiBI theory 

The linear dynamics of EiBI theory is equivalent to 
linearized Einstein's theory. To prove this, we expand 
the field equations (jU)-© around a vacuum Minkowski 
background, 

q a b = Vab + 5q ab , g a b = Vab + Sg ab , (15) 
q ab = if b - Sq ab , g ab = if b - Sg ab , (16) 
T a b = ST ab , T ab = ST ab , (17) 

where indeces are raised and lowered by the Minkowski 
metric, rj ab - In order for the Minkowski metric to be solu- 
tion, we set A = 1 in the field equations. Linearization of 
Eq. (j4|) follows immcdiatly from the standard linearized 
Ricci tensor, 

Sq a b-Sg ab = n5R ab = — {5q c a , bc + ^b.ac - OSq a b - 5q a b) , 

(18) 

where a coma denotes partial derivative with respect to 
the flat background and Sq is the trace of 5q a b- In order 
to linearize Eq. ([5]), we recall that, at linear order, |— g\ = 
1 + 5g, so we get 



jib 



5q ab -^{Sq- Sg) = 5g ab + nST ab 



(19) 



By taking the trace of the equation above, we get Sq 
Sg = —k5T, which can be substituted back into Eq. (fT9 
Finally Eq. (|18p can be written as 



5R ab = 5T ab - -^-ST , 



(20) 



which does not depend on k. This equation has exactly 
the same form as in GR, but the Ricci tensor is written 
in terms of q ab . However, at linear order, we note that 



8T n , 



5T n 



(5* 



M 



5^ 



ST ab 



M 



Sg 



Im i 



(21) 



where the second term vanishes when evaluated in vac- 
uum, since matter fields are coupled to the metric. 
Therefore, ST ab and ST do not depend explicitly on the 
metric and Eq. ([20]) is exactly equivalent to the linearized 
Einstein equations to all orders in k, but for the auxiliary 
metric q a b- 

In particular, the linearized version of EiBI theory ad- 
mits the same differential structure as linearized GR. For 
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example the principal symbols, i.e. terms in the lin- 
ear equations involving highest order derivatives, are the 
same as GR. Although a detailed analysis of the well- 
posedness of the theory is beyond our scope, these re- 
sults suggest that there exists a well-posed formulation 
of EiBI theory, similarly to GR. For a unitarity analysis 
of general BI gravity theories see Ref. [4(| . 

Finally, in the small K limit, our initial assumption of 
Minkowski background can be relaxed. In this case, sim- 
ilar conclusion about the linear structure of EiBI around 
any background metric can be drawn. 



III. NON-RELATIVISTIC STELLAR COLLAPSE 



where u is the fluid velocity, p is the density and m is the 
mass function. Note that the second equation above is 
vectorial, Q is the viscosity tensor and V • Q is a vector. 
The viscosity tensor is defined as [44| 



Q = £ 2 pV ■ u 



Vu - -V • u 



(25) 



if V • u < 0, otherwise Q = 0. In the equation above 
Vu = (djUi + diUj)/2 and e is the unit tensor. Note that 
this equation only contains scalar invariant quantities so 
that can be directly specialized to the spherically sym- 
metric case, contrarily to Eq. (8) in Ref. [43[, which is 
only valid in Cartesian coordinates. 
In indicial form, 



In this section, we describe the non-relativistic collapse 
of dust in EiBI theory. Hereafter, we focus on asymptot- 
ically flat solutions, setting A = 1. The collapse of in- 
coherent dust in the Newtonian limit shares many prop- 
erties with its relativistic analogue [4J [42j]. Here, we 
shall solve the hydrodynamics equations in the case of 
spherical symmetry, showing that the end-state of the 
1 + 1 evolution is the pressureless star found in Ref. (26j 
and described in detail in Section IIVI below. This con- 
firms and extends the analytical argument presented in 
Ref. [26| . However, at the end of this section we show 
that our analytical computation in Ref. (26j is partially 
flawed, due to a typo in the series expansion close to the 
origin, and we shall discuss why the numerical results 
are still in agreement with the overall picture given in 
Ref. [H. 

The hydrodynamical equations can be equivalently 
solved using either the Eulerian or the Lagrangian ap- 
proach. In the next two sections, we shall briefly review 
these formulations, together with a standard procedure 
to avoid shock wave formation, and the modifications re- 
lated to the problem at hand. 



A. Eulerian formulation 

In the non-relativistic limit, the collapse of non- 
interacting (P = 0) particles is governed by the mass 
conservation (continuity equation) and momentum con- 
servation (Euler equation), the latter being modified in 
EiBI gravity due to Eq. (ft]). 

Following the seminal work [43|, we supplement the 
hydrodynamics equations by an artificial viscosity term, 
in order to avoid divergences due to shock wave formation 
during the evolution. In the Eulerian formulation, the 
relevant equations in the pressureless case read [U |44[ 



dp „ 

<9u _ „ x 
— + u • Vu + V$ 

at 



- pV • u 


= 0, 


(22) 


V Q 
h 

P 




(23) 


J d 3 xp 


= o, 


(24) 



Qij = ( 2 pV ■ u 



djUj + diUj 



—V • u 

3 



(26) 



where I is a constant with dimension of length. 

We are interested in the spherical symmetry case, 
where all the dynamical variable only depends on (t, r) 
and the only non- vanishing component of vectorial quan- 
tities above is on the radial direction, e.g. u = 
(u(t, r), 0, 0). In this case only the radial component of 
Eq. (f2"3")l is nontrivial and [V • Q] r = Vq, where qt = Qi r . 
Due to the symmetry, the vector q = (q(t, r), 0, 0), and 



q(t, r) = £ 2 p{t,r)V ■ u 



u'(t,r) 



1. 



(27) 



Finally, in spherical symmetry, V ■ u = u' + 2u/r and the 
artificial viscosity term in Eq. (f2"3")l is V • q = q' + 2q/r, 
so that we are left with the set of PDEs 



du 
~dt 



dp , 

, Gm 
uu + — — 



pu 



—pu 

r 



0, 



P 



m 



47rr 2 p 



0, 



(28) 

(29) 
(30) 



where a dot and a prime denote derivatives with respect 
to t and r, respectively. The viscosity term reads 



U 2 

Q = — [rp (u + 2ru') u" 
3r z 

-(u- ru') (3pu' + (2u + ru') p')} 



(31) 



if u' + 2u/r < 0, otherwise Q = 0. In practice, I = cAr, 
where Ar is the grid spacing in the radial direction and 
c is some dimensionless constant. This allows the shock 
to be "smeared" in a region of width £, while leaving the 
rest of the dynamics unchanged. In our simulations, we 
have checked that the results do not depend on the pre- 
cise value of c. Actually, the results do not even depend 
on the details of the artificial viscosity term, as long as 
it introduces a smearing effect on the shocks while keep- 
ing the rest of the physics unaltered. We have explicitly 
checked this fact, by comparing test simulations with dif- 
ferent artificial terms. 
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Finally, one can derive a Bernoulli-like relation by mul- 
tiplying Eq. (J231) by u(t,r), leading to 

dtir^ + V=m- (32) 

where we neglected the artificial viscosity term and where 
d/dt is the total time derivative. This expression is valid 
locally and the term on the right hand side takes into 
account the local variation due to matter flow. Indeed in 
the pure Newtonian case, this term can be rewritten as 
a surface integral of the density current. 



B. Lagrangian formulation 



C. Results 

We have solved the initial- value problem defined by 
the above system of PDEs using the method of lines, in 
which the radial dimension is discrctized and the result- 
ing system of coupled ODEs is integrated in time with 
standard methods. The initial static profiles we consid- 
ered are presented in Table HI 





p(0,r) 


u{0,r) 


Profile I 







Profile II 


2 

sin zur ^—-ar 

■UJT 






TABLE I. Initial density and velocity profiles of our simula- 
tions. 



The basic Lagrangian equations are given in [43[ . They 
are expressed in terms of the comoving volume V(t, x) = 
^VX(i, a;) and the velocity field u = dX g^ , where 
Po(x) is the initial density. The continuity equation for 
the volume is given by 

p (x)V(t,x)=X(t,x), (33) 

which, in terms of the density p(t,x) = 1/V(t,x), reads 



p = -p^r, 



u'(t,x) du(t,x) 



X'(t,x) P 'dX{t,x) 



(34) 



If X(t,x) = R(t,x) is a radial variable in spherical 
coordinates, this equation becomes 

2 , . pit, x)u'(t, x) 
R{t,x) R'(t,x) 

Note that the equation above can be simply obtained 
from the corresponding Eulerian formulation by replac- 
ing the Eulerian time derivative by the Lagrangian time 
derivative according to 



D d 



(36) 



where D / Dt is the time derivative in the Lagrangian for- 
mulation. 

Following the same principle, the equation for the ve- 
locity field is given by 

Gm(t,x) p'(t,x) F q (t,x) 

u = ; 7^ K — 77 r H , , , (37) 

R(t,x) 2 R'(t,x) p (x) ' V ; 

where F q (t, x) is the artificial viscosity term given by [43j | 
which explicitly reads 



p = d_ f p (x)(cAx) 2 u'(t,x)\u'(t,x)\ 
q dx\ R'(t,x) 

where Ax is the grid spacing and c is a constant. 



(38) 



Profile I is a classical exponential density profile, while 
Profile II is a deformation of the pressureless star found in 
Ref. [H] (cf. Eq. (@SJ below). For both class of profiles, 
a is a free parameter related to the slope of the initial 
density profile at the center. 

For testing purposes, we wrote two independent codes 
that solve the hydrodynamical equations in the Eulerian 
and in the Lagrangian formulation, respectively. Once 
convergence is reached, the results of the two codes agree 
very well. However, the simulations in the Lagrangian 
formulation generically need a lower resolution, due to 
the fact that the radial grid evolves along with the fluid 
collapse. Indeed, the Lagrangian formulation is analog 
to the comoving frame generically adopted in relativis- 
tic collapse simulations. Hence, we shall present results 
obtained using this formulation. Typically, we used a 
non-uniform grid with 4 x 10 3 points in the spacial direc- 
tion, which guarantees convergence of the results for all 
values of k taken into account. As expected, we find that, 
for a given value of the viscosity constant c, convergence 
and stability of the numerical results is reached provided 
the mesh is sufficiently fine. 

In standard Newtonian gravity, n = 0, the hydrody- 
namics equations can be solved analytically for a con- 
stant density profile and they correspond to the relativis- 
tic Oppenheimcr-Snyder solution [42j. In that case, the 
dust collapses in a finite time tc = n\/R 3 / (8M), where 
R and M are the initial radius and the total mass of the 
spherical dust configuration. This analytical results is 
also valid when non-constant initial density profiles are 
considered (26|. Our simulations reproduce this result 
with very good precision. Furthermore, for any negative 
value of k, we found a qualitatively similar behavior, i.e. 
the fluid collapses to a singular state in a finite time. For 
k < 0, the time of collapse decreases with |«|. 

On the other hand, the behavior for any k > is dras- 
tically different, as shown in Figs. Q] and [2j We present 
results for np c = 0.1,1 and obtained using Profile I in 
Table (|T|, but different choices of n > 0, different initial 
profiles or different values of a, would give qualitatively 
similar results. In Fig.Q]wc show the central density as a 
function of time for two different values of n. The central 
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FIG. 1. Central density as a function of time normalized by its initial value. Left panel: np c =0.1. Right panel: np c = 1. 
Dashed lines denotes the central density of the pressureless star with the same mass (cf. Sec. IIV K\ . The oscillation period 
agrees very well with the fundamental period of proper oscillation of these solutions (cf. Sec. IIV Dl and Fig. @) 




FIG. 2. Radial profiles of the density (left panel) and of the velocity (right panel) for different instants and for np c — 1. The 
dashed vertical line denotes the radius of the pressureless star (cf. Sec. (|IV A|l ) with the same mass and the dashed black thin 
line in the left panel denotes the pressureless star profile. 



density is always finite and oscillates around a constant 
value at late time. As we shall show in Sec. ([IV A[) . the 
mean value of the oscillations (denoted by an horizon- 
tal dashed line in Fig. [T]) is precisely the central den- 
sity of a pressureless star (2(| (cf. Sec. (jIV A[0 with the 
same mass. In Fig. [2] we show the density and the ve- 
locity radial profiles at different instants. The would-be 
shocks in the velocity profile (the steep region shown in 
the right panel of Fig. [2]) propagates towards the exte- 
rior of the fluid without developing a discontinuity, which 
could not be resolved due to finite grid spacing. The arti- 
ficial term discussed above smears this discontinuity and 
ensures stability of the numerical simulations. 

Our simulations suggest that, even for an arbitrarily 
small value of k, this oscillatory behavior would continue 
indefinitely. This is perfectly consistent with the fact 
that the hydrodynamical equations do not include any 
dissipative term, so that energy is conserved during the 
evolution. Nevertheless, our results give strong evidences 
that, for generic (static) initial profiles, the end-point of 
the dust collapse is precisely the pressureless star. In- 
deed, we expect that any dissipative term, which has 
to be included in realistic situations, would quench the 



oscillations and drift the system toward a static config- 
uration, which is precisely the pressureless star. This 
picture is also confirmed by Fig. [3j which shows the pe- 
riod of oscillation, as a function of k, compared with the 
fundamental period of proper oscillation of a pressureless 
star (cf. Scc. lIVDl) . 

Finally, as shown in Eq. HU, the pressureless collapse 
in the modified non-rclativistic theory is equivalent to 
the collapse of a polytropic fluid with P(p) = up 2 / 8 in 
Newtonian theory. Hence, our results are consistent with 
the fact that, in the latter case, the final state is a regular 
polytropic star with polytropic index n = 1 (see e.g. [45j 
for more realistic simulations of star formation in New- 
tonian gravity). 

1. On the analytical method of Ref. fizdl l 

In Ref. [1^ , we developed an approximated method to 
solve the collapse equation analytically for any n. Here, 
we point out a typo in that computation and we discuss 
why this does not affect the final result, as we have proved 
in the previous section by solving the collapse equation 
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FIG. 3. Period of the oscillations around the pressureless 
static configuration as a function of n and with gaussian initial 
profile (cr = 10). The straight line is the prediction from the 
linear stability analysis (cf. Eq. (|61|> ) whereas the markers 
are extracted from the numerical simulations. 



numerically. 

Let us review the procedure adopted in Ref. [26| . Due 
to the spherical symmetry, we expect the singularity to 
initially form at the center r = 0. The collapse equations 
are then solved by expanding the dynamical variables 
close to the center, 



P 
u 



Po(t) 
u (t) 



Pi(t)r + p 2 (t)r 2 - 
-u x {t)r + u 2 (t)r 2 



0(r 3 ) , 
-u 3 (t)r 3 



0(r A 



and solving a system of ODEs for pi{t) and Ui{i). At first 
order, one finds pi{t) = uo(t) = u 2 (t) = and 



«i(*) = - 



Po(t) 
3po W 



At second order, we get 



Po(t) 



3 Po (t) 2 

hit) 



3 

~ 2 
5 Po (t) 



Kp 2 {t)-8wGp {t)=0 
Po(t) 



P2 (t) 3/>o(t) P2(t) 



u 3 (t) =0. 



(39) 

(40) 
(41) 



The results we presented in Ref. [261 ] originated by er- 
roneously omitting the term proportional to u 3 (t) in 

Eq. gTJ. Indeed, if u 3 (t) = 0, then p 2 {t) oc p 5 /3 (t), 
which can be substituted in Eq. (|40[) to obtain a single 
non-linear ODE for po(t). The solution of that equation 
is given in the equation below Eq. (7) of Ref. [26j . 

However, in general u 3 (t) is non- vanishing and 
Eqs. (|40[) and (|41j) are not sufficient to solve for the three 
variables 7/3, po and p 2 - In fact, it can be easily proved 
that, at any order in the series expansion, the number of 
unknown functions is always larger than the number of 
differential equations. This is precisely due to the extra 
derivative term in the Poisson equation ([1]). Indeed, in 
the standard case when k = 0, Eq. (|40|) decouples and 
can be solved in the usual way. As a matter of fact, when 
k ^ 0, the collapse equations cannot be solved exactly by 
this simple method, as erroneously stated in Ref. |26| . 



Nevertheless, as we have showed above, our fully nu- 
merical simulations not only confirm the results previ- 
ously reported, i.e. that the collapse does not lead to 
any singularity, but also that the final state is the pres- 
sureless solution reported in Ref. J2f|. Furthermore, the 
oscillatory behavior found in the numerical simulation 
qualitatively match the oscillatory behavior presented in 
Fig. 1 of Ref. H|. 

Thus, it is relevant to understand whether the analyt- 
ical method present above, although flawed, can actually 
reproduce the full solution in some limit. This is cer- 
tainly true when po(t)u 3 (t) -C p 2 (t) at any time. The 
latter condition is indeed satisfied when the velocity and 
density profiles are close enough to a late-time evolving 
configuration, for example when the configuration oscil- 
lates around the static configuration. Indeed in this case, 
Po (t) = p s + 0(e), p 2 = p s 2 + 0(e) while u 3 (t) = 0(e), 
where p$, p| are the constant central density and second 
derivative at the center of the static pressureless config- 
uration and e is a small number. It follows that in this 
case, the combination pou 3 is one order smaller than p 2 . 



IV. STARS IN THE NON-RELATIVISTIC LIMIT 

In this section, we discuss non-relativistic stellar mod- 
els in the theory defined by © . Our results apply to any 
relativistic theory which reduces to Eq. ([T]) in the non- 
relativistic regime. We shall solve the hydrostatic equi- 
librium equation, supplemented by the standard mass 
conservation, dm/dr = Airr 2 p(r), and an EOS. Requiring 
spherical symmetry, the hydrostatic equilibrium equation 
reads (Hi 



dP _ Gm(r)p 
dr 



(42) 



Clearly, equilibrium configurations in EiBI gravity are 
different from the standard Newtonian ones in presence 
of non-trivial density profiles. Note that the equation 
above can be written in a more evocative form as 



dP m(r)p(r) 



(43) 



where we have defined an "effective" Newton's constant 

(44) 



Ge ff (r)=G + r V(r) 



4 m(r) 

Since p'(r) < inside a star, G e g ^ G when n 0. In 
Ref. [331 ] . the modified hydrostatic equilibrium equation 
has been used to construct realistic solar models and put 
the first observational constraint on k (see also the recent 
Ref. (35| for more stringent constraints). 

Here, we shall mainly focus on polytropic models with 
EOS of the form P = P(p). Using the mass conservation 
equation, we can rewrite Eq. (|42p as 



1 



— P 

P 



jdP(p) 
dp 



4nGp+- \p" + 2t- ) =0, (45) 
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where a prime denotes derivative with respect to r. Once 
the EOS is fixed, Eq. (|4"5j) is a second order ODE for p(r), 
which can be solved by imposing regularity condition at 
the center, i.e. 



p(r) 



n=0 



0. 



(46) 



In the equation above, pi are constants which can all 
be expressed in terms of the only free parameter, po, by 
solving Eq. ([35]) pcrturbatively close to the center. For 
example, 



Pi 



SGirpl 



3(4F'(p ) + K P o) 



(47) 



whereas pi = for any odd /. It is easy to see that stellar 
models only exist when p 2 < 0, which implies a lower 
bound, np > —AP'(p). As we shall see, this limit holds 
also qualitativel y in the relativistic limit and in more 
realistic models 1331] . 



A. Non- relativistic pressureless stars 

In classical Newtonian gravity, a gas of non-interacting 
particle cannot support self-graviting configurations and 
it will inevitably collapse due to its self-gravity attrac- 
tion. However, in theories in which the Poisson equation 
reads as in Eq. ([I]), a straightforward solution of the hy- 
drostatic equilibrium equation with P = and n > 
reads (H 



p{r) 



sm(zur) 



w = 4 



The radius and mass of the star read 

4tt 2 p c 



7T 
W 



M 



(48) 



(49) 



respectively, so that for a given k > 0, the solution above 
is uniquely characterized by the central density or, cquiv- 
alcntly, by the mass. In particular, the radius of the 
compact object is uniquely determined by k. 

The existence of such solutions is remarkable not only 
because they are interesting models for self-gravitating 
dark matter Q, but also because they provide a non- 
trivial static and spherically symmetric solution of the 
hydrodynamics equations governing the collapse of dust. 
Indeed, as we discussed in Sect. IIII1 our numerical sim- 
ulation give strong evidences that, as an outcome of the 
collapse, any initial distribution of matter will eventu- 
ally relax towards the pressureless solution described by 
Eq. gSJ). 

Clearly, the generality of this result is only due to the 
fact that we considered the collapse of pressureless mat- 
ter, but it also suggests that, when taking into account 
realistic models in which the matter is described by a spe- 
cific EOS, the end-state of the collapse will be a modified 



compact star, rather than a singular state. A detailed 
analysis in that direction would be certainly interesting. 
For the time being, it is important to notice that these so- 
lutions are also linearly stable, as first shown in Ref. (26j 
and detailed in Sec. IIVDI 

We note here that the pressureless star (|48p is a sort 
of solitonic solution, since it solves the Poisson equation 
for any harmonic function $, due to the identically van- 
ishing of the right hand side of Eq. ([1]), i.e. G c ff = 0. In 
the interior, the Newtonian potential is constant and it 
matches continuously the vacuum potential M/r at the 
radius. Indeed, Eq. (QJ can be thought as a forced oscil- 
lator equation for p, where k would play the role of the 
inverse spring constant and $ is an external force. The 
pressureless star qualitatively corresponds to the solution 
of the oscillator with a constant external force. 



B. Newtonian polytropic models 

For a generic polytropic index n, the field equation 
must be solved numerically, imposing p ~ po+p 2 r 2 at the 
center. In this case stellar solutions only exist provided 
the following condition is satisfied 



K > -\k c \ = -4A'(1 + l/n)p\ 



-l + l/n 



(50) 



For k > 0, condition ([50|) is always fulfilled. In some 
cases the Lane-Emden equation obtained from Eq. (|42p 
can be solved analytically (46[. For instance if n = 1, 
P{p) = Kp 2 , and the solution reads as in (|48p . but with 
w = AyJ Gir/ (8K + k), so that it exists for k > — 8K 
and reduces to the pressureless case for K — 0. This is 
consistent with the effective pressure term defined in (|14l) . 
Indeed, the modified pressureless star solution (|48p with 
k = 8K is exactly the same as the n — 1 polytropic 
solution with n = 0. 



C. Modified Chandrasekhar model for white dwarfs 

As a relevant example of non-relativistic polytropic 
star, let us consider zero temperature non-rotating white 
dwarfs. The interior of a zero-temperature white dwarf 
is described by a relativistic EOS, P(p), which is para- 
metrically defined by [13] 



P(x) 



irmfc 5 p e mp 



, . Simile 3 u e mp ^ 

P( X ) = ^3 ^ 



x{2x 2 -3)v / T+a^ + 3sinh" 1 x 



where c and h are the speed of light and the Planck con- 
stant, m e and mp are the electron and the proton mass, 
respectively and p e is the molecular weight. In Fig. [4] 
we show the mass-radius relation for different values of 
k > obtained by integrating numerically the hydro- 
static equations. 
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FIG. 4. Mass-radius relation for zero temperature white 
dwarfs for different values of k (in units of a typical white 
dwarf density, pwD = 10 9 kg/m 3 ) in EiBI theory and com- 
pared to the standard result for k — 0. Notice that, for any 
K > 0, there is no maximum mass, but white dwarf models 
have a minimum radius given by Eq. (|55[) . 



When k = we recover the famous Chandrasckhar 
result: zero temperature white dwarfs have a maximum 
mass, M as 1.4M Q . However, for any n > 0, the be- 
havior is drastically different. White dwarf models can 
have an arbitrarily large mass, but instead there exists a 
minimum radius, which approximately scales as k 1 / 2 . 

We can use Landau's original argument to understand 
these results. Following Shapiro and Tcukolsky's exposi- 
tion (48|, let us consider a star of radius R composed of 
N fermions, each of mass m^. In the relativistic regime, 
the Fermi energy of the system reads 



E F = 



HcN 1 / 3 
~R 



(52) 



In EiBI theory, the gravitational energy per fermion is 
approximately 



E G 



GNm 2 b ZnNml 



R 



16ttP 3 



(53) 



where we approximated the star's density as p 
Nm,b/ {AttR? /3). Thus, the star's total energy read 



E — Ep + Eq = 



HcN 1 / 3 - GNml 'dnNml 



R 



WttR 3 



(54) 



For small N, the total energy is positive, and we can 
decrease it by increasing R. If moreover 



R> R, 



crit 



3« 
16ttG 



(55) 



then the total energy can become negative. At some 
point the electrons become non-relativistic and the 
purely Newtonian term dominates. It is easy to see that 
when this happens the energy tends to zero when R goes 
to infinity. Thus, there is a point of minimum energy, 
and it is delimited by the R = P C rit curve (|55|) . Note 



that the scaling i? cr it ~ \[k is consistent with the recent 
computation of the Jeans length in Ref. [35|. 

On the other hand, even if E < initially, for suffi- 
ciently small R condition (|55p will be violated: the EiBI 
term dominates at sufficiently small radii, and collapse is 
never energetically favored. Relativistic stellar models of 
white dwarfs are described in Sec. IV Dl 



D. Stability in the non-relativistic limit 

We now discuss stability of the Newtonian configura- 
tions against radial perturbations [48j , expanding the dis- 
cussion in Ref. [26|]. The equations governing the stellar 
dynamics are the Poisson equation ([T]) together with the 
continuity equation and the momentum equation. These 
are given in Eqs. (|22|) - ([24|) by replacing the artificial vis- 
cosity term by a physical pressure term, i.e. V • Q — > VP. 
In those equations, the hydrostatic equilibrium is recov- 
ered when u = 0. 

We shall focus on spherically symmetric models. The 
perturbed Euler equation reads 



A 



du 
~dl 



EL 

p 



= 0. 



(56) 



where d/dt = d/dt + u ■ V is the total time derivative, 
A = S + £ • V and £ is the Lagrangian displacement. 
In order to compute this equation explicitly we use the 
following dH 



Ap 

Sp 

AP 



-4w, 

-47rp£+ j{5p)' = -47rp£ 

P A P 
Pi , 

P 



(r 2 pZ)' 



where the last equation above defines the adiabatic index 
of the perturbations, 7. It is then straightforward to 
obtain the following 



4 'f 

P 



d 2l 
dt 2 ' 

Ap D , , A(P') 



= 2 ±P> 



(AP)' 



rp 
1 



rp p 
A($') = (A$)' - = + £$' 

= (5$)' + £V 2 $- 



2C 
pr 



P' 



<p -£,p 



where we have used Ad/dr = d/drA — ^'d/dr and 
the background equations. Finally, using the equations 
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above, Eq. (|56|) reads 

-~ip 



P 

K 

+ 4 



pr 



0. 



(57) 



Assuming a time dependence e lw *, the modified eigen- 
value equation reads (2rjj 



4gP' «/> 
r 4 



(58) 

This equation must be solved for £. requiring the follow- 
ing boundary conditions [48| 



£(0) = 0, AP(i?) = - 7 -(r 2 O' = 0, 



(59) 



the latter being equivalent to requiring regularity of £ at 
the radius. An instability corresponds to an eigenmodc 
with lj 2 < 0. 



Perturbations of pressureless stars 



For P = and p given by Eq. (|48|). the eigenvalue 
equation (|58[) simplifies. In particular it does not depend 
on the index of perturbations 7, but only on k. Thus, for 
a given k, there is one fundamental mode. By integrating 
Eq. (|58p numerically and imposing Eqs. (|59|) . we find 



ap, 



1/2 



(60) 



where a ~ 2.1866, independently from k. 

This gives the characteristic period of oscillation 



AT 
~W 



7T 5 / 4 / K \ 3 / 4 

~2^~ KM 2 ) 



(61) 



which is shown in Fig. [3] and compared with the period 
of oscillation of our numerical solutions around a pres- 
sureless star with same mass. The results of the linear 
analysis perfectly agrees with the simulations, giving fur- 
ther support for the end-state of the collapse. 

In addition, we found no unstable modes of the pres- 
sureless star, confirming that in the modified Newtonian 
theory these object are linearly stable (26| . 



2. Newtonian polytropic stars 

In Newtonian gravity, polytropic models with 7 = 4/3 
are marginally stable for any polytropic index n [48J . In 
our case, these models are stable if n > and unstable 
if k < 0. A representative example is shown in Fig. [5l 
where we show the fundamental mode as a function of 
k for the polytropic model with n = 1 and 7 = 4/3. 
For generic values of 7, positive values of k contribute to 
stabilize the models, while negative values work in the 
opposite direction. 
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FIG. 5. Fundamental mode of the polytropic star with n = 1 
and adiabatic index of perturbations 7 = 4/3. In the standard 
theory (« = 0) the star is marginally stable, whereas it is 
stable for k > and unstable for k < 0. 



V. RELATIVISTIC COMPACT STARS 

In this section, we study compact objects in the fully 
relativistic theory defined by the action ([2]). Let us now 
consider static and spherically symmetric perfect fluid 
stars. The metric ansatzen read 



q a bdx a dx b 
g ab dx a dx b 



-p(r)dt 2 + h(r)dr 2 + r 2 dtt 
-F(r)dt 2 +B(r)dr 2 + A(r) 



r 2 dn 2 



(62) 
.(63) 



where we have used the gauge freedom to fix the function 
in front of the spherical part of the auxiliary metric q. 

Notice that the field equations ([5]) are simply algebraic 
equations relating q and g. Inserting the ansatzen above 
into Eqs. ([5]), and solving for the coefficients of g in terms 
of q, leads to 



F = p(r)[l + Kp] 2 A 3 (r), 
B = h(r)A(r), 
A=[{l- K P){l + K p)}- 1/2 



(64) 
(65) 
(66) 



which correctly imply g^ v = q^ in vacuum. The dynam- 
ical field equations @ read 

k (p' (rph' + h (-4p + rgQ) - 2rhpp") 
F = p 7 - rT -_ , (67) 



Arh 2 p 



B = 



rFh 2 - K {ph' + hp') 
rhp 



K ( 2rti 4 2rp' 
4r z \ h h hp 



(68) 
(69) 



where F, B and A are given in Eqs. (|64p - (j66|) . The equa- 
tions above, supplied by an EOS in the form P = P(p), 
can be solved for p, h and P. Since matter is covariantly 
coupled with the metric g, the standard conservation of 
the stress-energy tensor follows, W a T ab = 0, where we 
recall that the covariant derivative is defined in terms of 
the physical metric, g a b- We consider perfect-fluid stars 
with energy density p(r) and pressure P(r) such that 



T 



ab 



rpab 
perfect fluid 



[p + P] u a u b + g ab P 



(70) 
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where the fluid four-velocity u a — (1/^/^,0,0,0). Thus, 
rather than using Eq. ()69|) . it is more convenient to use 
the conservation of the stress-energy tensor, 



2FP' + F'(P + p) = 0. 



(71) 



Notice that, although the equation above simply reads as 
in GR, F 1 introduces terms proportional to p 1 through 
Eq. (|64|) . Finally, in the limit k 0, the field equations 
reduce to that of GR with 8irG = 1. 



A. Integration 

We have integrated Eqs. (|67|) . (|68|) and (fTTj) imposing 
regularity conditions at the center of the star, where the 
following expansions hold 



p(r) r 
P(r) 



Po - 

- P 



p 2 r , 

f P-2T 2 



h(r)r. 
p(r) 



1 + h 2 r 2 , 

J Pc + P2T 2 



We can set po = 1 by a time reparametrization. Fur- 
thermore, using the field equations, the coefficients 
(P2, ^-2, P2) can be written in terms of P c and p c as follows 

2(1 - P c k){P c + p c )(l + np c ) ((1 + Atpe) 2 ^ - 1) 
2 - 3kA 

A = (P c k - 3k Pc - 4)(1 + KPc) - k(1 - kP c )(P c + p c )^ c , 



= 



(1 + K Pc )*Al 1 
3k 



/l 2 = 



2 + (1 + K Pe ) 2 ,4g - 3A e 
6k 



with p^, = dp(P c )/dP, A c = A(0) and where we have 
used p 2 = Pip'c- 

The series expansion of the field equations at the 
center of the star contains terms of the form A c = 
y/JT— kP c )(1 + Kp c ). Assuming p Cl P c > 0, k must sat- 
isfy two conditions in order to allow for self-gravitating 
objects: 



P c k < 1 , 
Pc\k\ < 1 



for k > , 
for k < . 



(72) 
(73) 



Hence, assuming NSs may reach a central density p c ~ 
10 18 kg m~ 3 and P c ~ 10 34 N m -2 , the bounds above 
strongly constrain the theory, \k\ < 1 m 5 kg _1 s^ 2 (see 
also the recent Ref. [35|]). Furthermore, it is easy to prove 
that compact objects only exist if Pi < 0, which requires 
kA < 0. This gives a further constraint depending on p c , 
P c and p' c . The form of the constraint is cumbersome, but 
it is qualitatively equivalent to Eq. (|50|) . In particular, 
the condition is always satisfied for k > 0. 



B. Israel-Darmois matching conditions 

The field equations are integrated outward up to the 
radius R, defined by the condition P(R) = 0, where 
we require the numerical solution to match the exact, 
and unique, vacuum Schwarzschild solution, F(r) = 



S(r)" 1 = p{r) = h{r)~ l = 1 - 2M/r, where M is the 
mass of the star. 

In GR, the mass is computed by imposing Darmois- 
Israel matching conditions [49[ at the radius, i.e. [gij] = 
and [Kij] = 0, where [...] is the jump across the surface, 
Kij is the extrinsic curvature tensor, and i,j = 0,2,3. 
These matching conditions comes from the Einstein equa- 
tions and the requirement of a well-defined 3— geometry. 
For a spherically symmetric spacetimes, this leads to 
[gtt] = = [g' tt ] and [g rr ] = = [g' rr ] through the field 
equations. 

The same procedure can be applied to EiBI theory, 
as we explain in detail in Appendix [Al In this case, the 
junction conditions read 



[<fe]=0, [ qij } = 0, [K l3 (q)} 



r Qij 



= 0, (74) 



where Kij(q) is the extrinsic curvature tensor for a hy- 
persurface of a manifold (Ai,q). 
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FIG. 6. Mass-radius relation for a relativistic pressureless 



stars. Results are normalized by po = 8 • 10 



m , which 



is a typical central density for NSs. The shaded region is 
excluded by causality, R > 2.9GM [5(| (cf. also the discussion 
at the end of Section |VE|). 



C. Relativistic pressureless stars 

The existence of Newtonian pressureless stars makes it 
relevant to investigate the existence of similar solutions 
in the full theory. To this purpose, we set P = 0. The 
conservation of the stress-energy tensor simply implies 
F(r) =const. For a given value of k, the solutions of the 
field equations then depend only on the central density. 

As shown in Fig. [6j for any value of k > 0, there exists 
a regular solution. These stars have a maximum com- 
pactness of GM/R ~ 0.3, which is roughly independent 
from k, a maximum masfp and a maximum radius which 
linearly increase with k. Furthermore, they have a posi- 
tive binding energy, fh/AI — 1, where fh is the baryonic 



1 We thank Jan Stcinhoff for having pointed out this property. 
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mass, corresponding to the energy that the system would 
have if all baryons were dispersed to infinity. 

Of course, such compact objects do not exist in GR, 
while they exist in EiBI theory because k > introduces 
a repulsive gravity contribution. Interestingly, the EOS 
for dark matter particles is approximately P = 0. Hence, 
in this theory self-gravitating objects, purely made by 
dark matter, can exist and may reach the typical com- 
pactness and mass of most compact NSs. 



D. Relativistic zero-temperature white dwarfs 

With the modified Tolman-Oppcnheimer-Volkoff equa- 
tion at hand, we can also study relativistic models of zero- 
temperature white dwarfs, described by the EOS ([51~j) . 
Note that this EOS can be explicitly written as P = 
P(p), whereas the inverted function p = p(P) must be 
written in parametric form. However, we have already 
discussed that in EiBI gravity the field equation (fTTj) also 
contains p' terms, so that it can be used equivalently as 
an evolution equation for p or for P. The former choice is 
more convenient, since we can substitute P(p) explicitly 
in the field equations and solve them for p and for the 
metric fields. 

Some results are shown in Fig. [7J We focus on k > 0, 
since this is the most interesting region of the parameter 
space. The left panel of Fig. [7J shows the mass-radius re- 
lation for a zero-temperature white dwarf in the full EiBI 
theory. Note that, for any k > 0, there exists a maximum 
mass. In the n — >• limit, the classical Chandrasekhar 
result, M max ~ 1.4M©, is recovered. This is shown in 
the left panel of Fig. [TJby an dotted horizontal line. 

When k is non-vanishing, the maximum mass can be 
much larger and it can easily exceed 10 2 M Q . Therefore, 
the absence of white dwarfs with mass M ^> Mq can 
be used to constrain the coupling k. Curves in Fig. [7] 
terminate when the central density is so high that condi- 
tion (|72p is not fulfilled, but the maximum mass gener- 
ically exists for smaller central densities. Although the 
mass can be quite large, the compactness GM /R of the 
object is limited, as shown in the right panel of Fig. [7] 



for rii-i < n < m, respectively Here mjn is the rest- 
mass density and p is the total energy density. Given 
a set of Ti and an initial Pi = P(mtni), the values of 
Ki are obtained by imposing continuity. In Ref. (5lj , the 
best fits for Ti for a three division piecewise are provided 
for many tabulated EOS. Here, we shall focus on the FPS 
EOS, for which [U 

r x = 2.985, T 2 = 2.863, T 3 = 2.600 , (76) 

and Pi = 6.653 • 10 34 g/(cm s 2 ). Piecewise polytropes 
can approximated tabulated EOS extremely well. With 
the choice above, the fit to the tabulated FPS EOS gives 
a residual 0.0050 and an error for the maximum mass 
and momentum of inertia of order 0.01%. Some results 
obtained using this EOS are shown in Fig. [8j 

Notice that, within GR, a standard EOS like FPS 
would be ruled out by the recent observation of a NS 
with M = (1.97 ± 0.04)M o [H3| (denoted by an horizon- 
tal band in Fig. [8j In EiBI gravity the maximum mass 
of a NS can be much larger than in GR, thus such obser- 
vation can be accommodated without invoking a stiffcr 
EOS. Curiously, no constraint on k comes from causal- 
ity, R > 2.9GM, [5(|. The latter is always satisfied even 
for very large values of k. This is related to the exis- 
tence of a maximum compactness, GM/ R < 0.3, which 
is independent from the EOS (see e.g. Figs. [6] and [7|). 
This interesting property have a general interpretation 
in terms of effective stress-energy tensor [39| . 

F. Slowly rotating models 

Slowly rotating stars can be constructed from the cor- 
responding static solutions (53|. At first order in the 
rotation, g tLp = —((r)r 2 sin 2 9, qt v = —t](r)r 2 sin 2 and 
the stress-energy tensor for a rotating fluid can be built 
from from Eq. (|70p with 

tt° = {« t ,0,0 i nu t } , u*= yj~(g tt + 2n tip + n 2 gipip ), 

where H is the angular velocity of the fluid. 

Solving the field equation at order 77, fi), we find 
two equations for £ and r\ 



E. Realistic nuclear-physics based EOS 

Using tabulated, nuclear-based EOS to study NSs in 
EiBI theory is numerically challenging. The presence of 
derivatives of the matter fields requires interpolation of 
the EOS, which can be imprecise, specially close to the 
radius, where P ~ 0. To avoid this problem, we have 
used a piecewise polytropic EOS [5l| . which can be con- 
structed analytically as follows. For a set of dividing 
baryon number densities, uq < rii < 712- ■■, we define 



p(P) = m b n + 



1 



r<-i 



-p. 



P = Ki(m b n) 



(75) 



ri(r)= \Jj^ [(1 - «P)C(r) + ^A(P + p)} 

0= 4h 2 [r 2 C, + (-r 2 + k) 77] + rnh' (2rj + rrf) 

-nh r— (2rj — rrj') + 2 (2rj + 4rrj' + r 2 rj") 
P 



The first equation above is an algebraic relation be- 
tween 77 and £. Substituting this into the second equa- 
tion above gives a second order ODE for (, which has 
to be solved by imposing regularity of £ at the center, 
C(0) = C c , C'(0) = and requiring continuity at the stel- 
lar radius. At infinity, the asymptotic behavior of £ reads 



14 




FIG. 7. Zero-temperature relativistic white dwarfs in EiBI theory. Left panel: mass-radius relation. The horizonal line denotes 
Chandrasekhar limit, M ~ IAMq. Right panel: compactness as a function of the central density. The horizontal line denotes 
the maximum compactness allowed by causality, GM/R < 0.35 [H^l (cf. also the discussion at the end of Section fV Efl . Results 
are expressed in units of a typical white dwarf density, pwD = 10 9 kg/m 3 . Curves terminate when condition (|72p is not fulfilled. 




FIG. 8. Compact stars in EiBI theory with FPS EOS obtained by fitting a piecewise polytrope model for different values 
of /t. Left panel: mass as a function of the central baryonic density pb- Right panel: mass-radius relation. Inset: binding 
energy as a function of pb. Results are normalized by po — 8 ■ 10 17 kg m~ 3 , which is a typical central density for NSs. Curves 
terminate when conditions (|72|l or (|73[) are not fulfilled. The horizontal band denotes the recent observation of a NS with 
M = (1.97 ± 0.04)M© [H, whereas the shaded region is excluded by causality, R > 2.9GM [H} 



C — ^ Coo + 2J/r 3 , where J is the angular momentum. For 
asymptotically flat solutions wc must impose Coo = 0. In 
Fig. Owe show the moment of inertia / = J /CI as a func- 
tion of the stellar mass for slowly rotating stellar models 
obtained with the piecewise FPS EOS discussed above. 

Note that Figs. El and p are qualitatively very similar 
to Figs. 3 and 4 in Ref. [2g|, which have been obtained 
using a simple polytropic EOS. 



VI. CONCLUSION 

EiBI theory is a viable modified gravity which intro- 
duces a highly non-trivial coupling to matter, while be- 
ing equivalent to Einstein's gravity in vacuum. In this 
theory, the metric tensor is minimally coupled to matter, 
thus enforcing the usual stress-energy tensor conservation 



and geodesic motion. Nevertheless, the non-local Born- 
Infeld structure of the theory affects the matter sector by 
introducing higher order corrections in the matter fields 
and in their gradients. 

Even at non-relativistic level, these extra terms allow 
for a new interesting phenomenology, which is compat- 
ible with current experiments. We have shown that in 
a class of non-relativistic theories - including the non- 
relativistic limit of EiBI gravity - described by a mod- 
ified Poisson equation ([T]) with k > 0, the collapse of 
non-interacting particles produces a regular, pressureless 
stars, rather than a singular state. We have discussed the 
stability of these and other non-relativistic stellar models, 
showing that pressureless stars are stable against linear 
radial perturbations and that positive values of k gcner- 
ically tend to stabilize the star. We have also extended 
original Chandrasekhar's analysis, showing that, for an 
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FIG. 9. Moment of inertia for FPS EOS obtained with a 
piecewise polytropic models as a function of the stellar mass 
for different values of k. Results are normalized by po — 
8 ■ 10 17 kg m -3 , which is a typical central density for NSs. 



arbitrarily small value of K > 0, zero-temperature white 
dwarfs and NSs do not have a maximum mass, but pos- 
sess a minimum radius R ~ k 1 / 2 . 

In the fully relativistic EiBI theory, we have con- 
structed several static and slowly rotating compact ob- 
jects: relativist pressureless stars which do not exist in 
GR, zero-temperature white dwarfs with relativistic EOS 
and nuclear-physics base models of NSs. All these solu- 
tions have a maximum compactness (always smaller than 
the causality bound) and a maximum mass, which generi- 
cally increases as a function of k and can easily be 10— 10 3 
larger than in GR. This suggests that, in the fully rela- 
tivistic theory, the collapse of very massive stars may still 
form black holes, but this may require collapsing stars of 
very high mass, depending on the value of k. 

For these reasons, the most natural and important ex- 
tension of our work would be to study the relativistic 
stellar collapse and identify the final state. This would 
be particularly relevant to understand whether singular- 
ities may still be formed during the matter collapse in 
EiBI gravity, or if their formation is prevented, similarly 
to what happens in EiBI cosmology |25j . 

Near-future experiments (e.g. the proposed NICER 
mission) will investigate the interior of NSs, allowing for 
null tests of GR inside matter. Measurements of the NS 
mass, radius and moment of inertia will also allow to 



perform tests of alternative theories [54| such as EiBI and 
to constrain the coupling parameters of the theory j3fjj . 
Since NSs are the most extreme matter configurations in 
the universe, such tests would constrain the matter sector 
of gravitational theories to unprecedented levels. 
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Appendix A: Matching conditions 

In this appendix, we derive the matching condi- 
tions (|74[) . focusing on the spherically symmetric case. 
In GR, one imposes Darmois-Israel matching condi- 
tions HH, i.e. 



[<?«] = 0, [Ka] 



d r gij 



0. 



(Al) 



across any hypersurface £ of a manifold (M,g). In the 
formula above, [...] is the jump across the surface, Kij 
is the extrinsic curvature tensor for £ and i,j = 0,2,3. 
The matching conditions comes from the Einstein equa- 
tions and the requirement that £ has a well-defined 
3— geometry. For a spherically symmetric spacetimes, 
starting with a metric ansatz 



ds 2 = -F(r)dt 2 + B{r)dr 2 



(A2) 



the Israel conditions implies [F] = = [F' /B]. The first 
condition can be fulfilled by a time reparametrization, 
whereas the second one defines the mass of the spacctime. 
Incidentally, these two conditions together with Ein- 
stein's equations imply that all the metric functions are 
C 1 . Indeed, one can equivalently impose [F] = [B] = 0. 

The same reasoning can be applied to EiBI theory. 
From the ansatzen (f63|) - (|62p . we now have five metric 
functions F, B,p, h, A and two matter fields, P and p. 
Note that, inside matter, q a i> is the metric which de- 
fines the metric connection T and, in turn, the covari- 
ant derivatives and the curvature tensor R a b- Thus, one 
would expect the mass to be directly related to the aux- 
iliary metric q a b, and not to g a b- The junction condi- 
tions (|74p follow from the fact that R a b appearing in 
Eq. (j4]) is defined in terms of the auxiliary metric q. Last 
of Eqs. (|74l) is a consequence of the field equations, while 
the first two conditions arc imposed to ensure that any 
hypersurface (either of the manifold (AA,g) or of (Ai,q)) 
have a well-defined 3— geometry. 

The matching conditions (f74")) can be derived as fol- 
lows [55| . Starting from Eq. ([5]), we perform the integral 
along the proper distance, n, measured perpendicularly 
through £ 



lim/ dnJ=q~q ab = lim / dn\J^ (g ab - nT ab 



(A3) 
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Requiring the absence of discontinuities or delta tunc- it follows [iQj (q)] = to guarantee regularity of q. 
tions in g (to ensure a well-defined 3— geometry) and pro- 
vided that T ab does not contain delta function terms, it 
follows that q ab is regular at the interface. Hence, from 
Eq. (|4]), we obtain lim c _>o j^^dnRa^T) = 0. Finally, 
from the decomposition of the four-dimensional Ricci ten- 
sor, r!u , into the Ricci tensor of the 3- metric, 



R%)=R%-2K u K]+K\K ij -^-^^, (A4) 
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